data <- read.delim("~/Desktop/webtraffic.txt")
col_total <- colSums(data)
Traffic <- matrix(col_total, nrow = 9, ncol = 9, byrow = TRUE)
Traffic[9,1] = 1000
Traffic
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 0 447 553 0 0 0 0 0 0
## [2,] 0 23 230 321 0 0 0 0 63
## [3,] 0 167 43 520 0 0 0 0 96
## [4,] 0 0 0 44 158 312 247 0 124
## [5,] 0 0 0 0 22 52 90 127 218
## [6,] 0 0 0 0 67 21 0 294 97
## [7,] 0 0 0 0 0 94 7 185 58
## [8,] 0 0 0 0 262 0 0 30 344
## [9,] 1000 0 0 0 0 0 0 0 0
Directed Graph
The Markov Chain is irreducible because all states communicate with each other. The Markov Chain is ergodic beacuse it is recurrent and aperiodic.
row_total <- rowSums(Traffic)
P = Traffic / row_total
P
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0 0.44700000 0.55300000 0.00000000 0.0000000 0.00000000 0.00000000
## [2,] 0 0.03610675 0.36106750 0.50392465 0.0000000 0.00000000 0.00000000
## [3,] 0 0.20217918 0.05205811 0.62953995 0.0000000 0.00000000 0.00000000
## [4,] 0 0.00000000 0.00000000 0.04971751 0.1785311 0.35254237 0.27909605
## [5,] 0 0.00000000 0.00000000 0.00000000 0.0432220 0.10216110 0.17681729
## [6,] 0 0.00000000 0.00000000 0.00000000 0.1398747 0.04384134 0.00000000
## [7,] 0 0.00000000 0.00000000 0.00000000 0.0000000 0.27325581 0.02034884
## [8,] 0 0.00000000 0.00000000 0.00000000 0.4119497 0.00000000 0.00000000
## [9,] 1 0.00000000 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000
## [,8] [,9]
## [1,] 0.00000000 0.0000000
## [2,] 0.00000000 0.0989011
## [3,] 0.00000000 0.1162228
## [4,] 0.00000000 0.1401130
## [5,] 0.24950884 0.4282908
## [6,] 0.61377871 0.2025052
## [7,] 0.53779070 0.1686047
## [8,] 0.04716981 0.5408805
## [9,] 0.00000000 0.0000000
a <- c(1, rep(0,8))
prob5 = a %*% P %*% P %*% P %*%P %*% P
prob5[5]
## [1] 0.1315178
The probability of a vistor being on Page 5 after 5 clicks is 0.1315178
Q <- t(P) - diag(9)
Q[9, ] <- rep(1,9)
rhs <- c(rep(0,8),1)
Pi <- solve (Q, rhs)
Pi
## [1] 0.15832806 0.10085497 0.13077897 0.14012033 0.08058898 0.07583914 0.05446485
## [8] 0.10069664 0.15832806
B <- P[1:8, 1:8]
Q <- diag(8) - B
rhs <- c(0.1,2,3,5,5,3,3,2)
m <- solve(Q,rhs)
m[1]
## [1] 14.563
The average time a vistor spends on the website is 14.563
\[ n \geq \frac {\frac{1}{\lambda^2}} {{{(10^{-3}})^2}*0.01} \] \[ n \geq\frac{10^7}{\lambda^2} \]
set.seed(400)
n1 <- 10000000/ (1^2)
n2 <- 10000000/ (2^2)
n3 <- 10000000/ (4^2)
x1 <- runif(n1,0,1)
x2 <- runif(n2,0,1)
x3 <- runif(n3,0,1)
y1 <- -log(x1)/1
y2 <- -log(x2)/2
y3 <- -log(x3)/4
i1 <- sum(sin(y1)) / n1
i2 <- sum(sin(y2)/2) / n2
i3 <- sum(sin(y3)/4) / n3
i1
## [1] 0.499993
i2
## [1] 0.199819
i3
## [1] 0.05892617
Metropolis-Hastings
q <- function(x, lambda){
return(lambda * exp(-lambda * x))
}
f <- function(x, k=2, theta = 2){
return(x^(k-1)* exp(-x/theta))
}
x_list = c(1, rep(0,14999))
for(t in 0:14999){
if (t == 0){
x = 1
}else {
x = -log(runif(1,0,1)) / x_list[t+1]
}
a = f(x) * q(x_list[t+1], x) / (f(x_list[t+1]) * q(x,x_list[t+1]))
u = runif(1,0,1)
if (u <= a){
x_list[t+2] = x
}else
x_list[t+2] = x_list[t+1]
}
hist(x_list[seq(0,15000,100)], freq = FALSE, breaks = 30, main = "Distribution of sampling", xlab = "Sample Value", col = "RED")
sample <- rgamma(n=100, shape=2, rate=2)
hist(sample, breaks = 20, main = "Gamma Distribution", xlab = "Value")
acf(x_list[seq(5000,15000,100)])
The sample generated in Problem 3b is random due to the acf, and the shape of the distribution is matched Gamma distribution